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ABSTRACT 

We describe results from large numbers of iV-body simulations containing from 250 to 1000 
stars each. The distribution of stellar masses is a power law, and the systems are isolated. 
While the collapse of the core exhibits the expected segregation of different masses, we 
find that the post-collapse evolution is, at a first approximation, homologous. This is 
quite surprising because there is no reason for supposing that mass segregation should not 
continue to have a substantial effect on the evolution of the cluster. In fact the spatial 
distribution of the mean stellar mass is nearly static throughout the post-collapse regime, 
except for the overall expansion of the systems, and this helps to explain why the post- 
collapse evolution is nearly self-similar. Self-similarity is also exhibited by the distribution 
of anisotropy and the profile of departures from equipartition, which show little change 
during the post-collapse phase. The departures from energy equipartition and isotropy 
are small in the core and increase with radius. During post-collapse evolution massive 
stars (mainly) are removed from the system by binary activity. This effect dominates the 
preferential escape of low-mass stars due to standard two-body relaxation processes. 
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1 INTRODUCTION 

This is the third in a series of papers in which we investigate the evolution of stellar 
systems using iV-body simulations. This has often been done before, but one of the main 
drawbacks of most such studies is the poor statistical quality of much of the data. Even in a 
system with 1000 stars, the sampling noise in the results becomes overwhelming whenever 
measurements are taken from a small subset of the stars. Therefore it is difficult to carry 
out detailed quantitative studies of many interesting phenomena, e.g. the evolution of 
the core, the spatial distribution of anisotropy, the escape rate, and so on. This problem 
becomes much more acute whenever stars of different masses are present, if it is desired 
to investigate mass segregation, departures from equipartition, the differential escape rate, 
etc. Indeed this paper is concerned with precisely this kind of issue. 

1 Now at the N. Copernicus Astronomical Centre, Bartycka 12, 00-716 Warsaw, Poland 
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The main technical improvement in the present research is an extremely simple one. 
Instead of presenting results from one simulation, or only a small number, we combine sta- 
tistically the results from several dozen simulations. The intention is to improve the signal- 
to-noise of our results to the point at which reliable quantitative data can be presented, 
and compared with the results of simplified models (e.g. gaseous models, or Fokker-Planck 
models.) An alternative method of improving the statistical quality would be to simulate 
a single, larger system. An advantage of that approach is that it works in the direction 
of greater realism (if the intention is to understand the dynamical evolution of globular 
star clusters) . But a much greater increase in computational effort is required for a given 
improvement in signal-to-noise, compared with our strategy of simulating large numbers of 
more modest simulations (see Giersz & Heggie 1994a). This strategy is especially effective 
now that parallel computers have become readily available. 

The present paper differs from the first two papers in this series (Giersz & Heggie 
1994a, 1994b, hereafter referred to as Papers I and II) by dealing with the case of unequal 
stellar masses. This greatly extends the range of parameters to be explored, but we have 
restricted ourselves to parameters comparable to those in the well known Fokker-Planck 
survey carried out by Chernoff & Weinberg (1990). (One of the motivations for our research 
was the desire to test their results using A-body methods, which make fewer assumptions.) 
On the other hand we still consider systems which are isolated (without tidal effects) and 
in which the stellar masses are constant, i.e. we do not model systems in which the stars 
lose mass through internal evolution. Those two processes are included in later sets of 
models which we discuss in two further papers. 

The plan of the paper is as follows. We begin with a discussion of the initial param- 
eters which we adopted, such as the slope and range of the stellar mass function. Then 
subsequent sections discuss, in turn, the spatial evolution of the systems, and in particular 
the segregation of masses. Next we discuss the distribution of velocities, and in particular 
the extent to which equipartition is obeyed in the systems at different radii. Then there 
follow sections on the evolution of the core and of the binaries which mostly form there. 
Towards the end of the paper we discuss global parameters, e.g. the evolution of the total 
mass and energy, with a detailed description of the differential rate of escape (i.e. its 
dependence on stellar mass). The last section discusses one or two points of interest and 
sums up. 

2 ORGANISATION OF THE CALCULATIONS 

As in earlier papers of this series, all models reported here started as Plummer models. 
The initial mass spectrum was chosen to be a power law of index 2.5, i.e. dN(m) oc 
m _2 ' 5 (im, in a range such that m max /m m i n = 37.5. This choice was made in order to 
facilitate comparison (in later papers in this series) with results of Chernoff & Weinberg 
(1990), who used a very similar spectrum. 

Though the spectrum of masses is continuous, for purposes of analysis the masses 
were grouped into 10 classes, along the lines of those used by Chernoff & Weinberg (see 
Table 1 in this paper). (In their Fokker-Planck calculations it is impractical to include a 
continuous spectrum of masses. Therefore a finite number of masses are used, and after 
considerable experimentation Chernoff & Weinberg found that use of 10 masses yielded 
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results quite similar to those using larger numbers of masses.) 

The calculations were carried out on the Edinburgh Computing Surface (ECS), which 
was a transputer array containing about 400 processors (cf. Paper I). The code used was 
NBODY5 (Aarseth 1985). Three sets of calculations were carried out, in which the number 
of stars was N = 250, 500 and 1000, respectively. In each set, the respective number of 
models was 56, 56 and 40. Virtually all runs were completed to time t = 400. 

The organisation of the output was handled somewhat differently than in earlier papers 
in this series. The reason for this was fairly directly related to the use of a spectrum of 
masses, divided into the 10 classes described above. Because of the steepness of the mass 
spectrum, the average number of stars in the class with heaviest stellar masses was less 
then one for N = 250. Since we wished to measure Lagrangian radii for each mass group 
in a uniform manner, and consistently with earlier papers where data for the inner 1% 
Lagrangian radii were discussed, it was necessary to combine data from all models at the 
same iV-body time and also to combine data from three mass bins corresponding to the 
most massive stars. (Even so, with only 56 models in each series, there are still only of 
order 150 stars in the heaviest mass group). Since the calculations produce data at different 
rates, however, this approach meant that it was necessary to store data on each star from 
each model at each output time. In order to conserve space this was done by coding data in 
such a way as to allow it to be reconstructed to adequate accuracy in subsequent analysis. 

The data stored (to limited accuracy) for each star was as follows: 

1. mass; 

2. distance from the density centre (defined as in Casertano & Hut 1985, with a modifica- 
tion for the presence of a mass spectrum, described below); and 

3. radial and tangential velocity (the radial direction being that from the density centre). 
In addition, the following details of all regularised binaries, and all unregularised binaries 
with binding energy > 0.1/cT (where the current mean kinetic energy of all single stars is 
3&T/2) were output: 

4. masses of components; 

5. separation of components, and distance of binary centre of mass from the density centre; 

6. internal energy of the binary and its eccentricity; 

7. radius of neighbour sphere and number of neighbours (defined as in NBODY5); 

8. names of binary components. 

The typical size of an output file from 40 runs at N = 1000 for 400 output times was 
of the order of 125Mb. Analysis of such a file proved to be manageable, if awkward, even 
on a diskless workstation requiring access to files over a fairly busy ethernet. 

In addition, each model outputted some preanalysed data, as described for equal-mass 
models in Paper I. For example, the values of ten Lagrangian radii (defined below) were 
output. As in our earlier work, the different sets of calculations did not all output every 
item in the list in the appendix to Paper I: the data selected for output was changed 
slightly as our interests and experience developed. 

The previous discussion refers to "Lagrangian radii" , a term which requires somewhat 
careful definition in the context of these multi-mass models. First, these radii are referred 
to the "density centre", which is defined as follows: for each star the density is defined 
to be proportional to M 5 r^ 3 , where r$ is the distance to the fifth nearest neighbour (cf. 
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Casertano & Hut 1985) and M 5 is the total mass of the five nearest neighbours. Then the 
density centre is defined as in Aarseth & Heggie (1992). Incidentally the density centres of 
different models were found to gradually move away from the origin (whose position and 
velocity coincide initially with those of the barycentre) in different directions. The root 
mean square value of each coordinate of the density centre (averaged over all systems in 
a given series) increases nearly linearly with time in the post-collapse regime, roughly as 
0.05t for 50 < t < 200 for N = 250. 

With the density centre thus defined, a Lagrangian radius is defined to be the radius 
of a sphere centred on the density centre and containing a fixed fraction of the total 
bound mass of the system. (When this does not correspond to a whole number of stars, 
interpolation of ordered radii was used as in Paper I.) Because of mass segregation and 
the presence of a spectrum of masses, it can happen that the innermost radius, which 
corresponds to a mass fraction of 1%, may contain less than one star in a single model. 
For this reason the Lagrangian radii discussed in this paper were actually determined 
differently. For each output time, the radii of all stars (measured from the density centre 
of the model to which it belongs) were aggregated from all models in a given series, and 
then Lagrangian radii were determined as above from this much larger set of stellar data. 
It was hoped that this would improve the determination of the inner Lagrangian radii, and 
avoid the difficulties associated with mass segregation, as described above. 

3 SPATIAL EVOLUTION 

3.1 The entire mass distribution 

The essential facts about the evolution of the mass distribution are apparent from 
Fig.l. The collapse of the core is much more rapid than in the equal mass cases considered 
in Papers I and II (cf. Table 2). Indeed the ratio of t co u (between the values for unequal 
and equal masses) is a decreasing function of A even for the largest value of A that we 
have studied. At first sight this appears to contradict the result that would be expected 
from the standard theory of relaxation, which is that this ratio should be independent of 
A. However, it was argued long ago by Henon (1975) that the coefficient 7, which appears 
in expressions for the relaxation time in the Coulomb logarith In 7 A, is generally smaller 
in systems with unequal masses than in those with equal masses, by as much as a factor of 
5. Indeed his theory yields the result 7 ~ 0.044 for the mass function used in our models. 

Now we turn to the numerical data (Table 2). If it is assumed that the collapse time 
(as a function of A for unequal masses) is proportional to the initial half-mass relaxation 
time, then we find that the results of Table 2 are roughly consistent if 7 ~ 0.015, which is 
a factor 7 smaller than the value found in Paper I for systems with equal masses. 

An independent method of determining 7 is comparison of t co u with the results of 
a Fokker-Planck computation with the same initial conditions and mass spectrum. This 
has been done using an isotropic code (Inagaki & Wiyanto 1984), from which the value of 
tcoii was obtained by finding the time at which the central potential reached the minimum 
attained in the A-body results (Fig. 12). Scaling to the values of t co u in Table 2 yielded 
values of 7 in the range 0.016 <^ 7 <^ 0.026. 

Whatever the correct value of 7, we have seen that there are three independent lines 
of evidence for a substantially smaller value than in the case of equal masses. This implies 
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that considerable care is required if the results of models with small A are to be scaled 
to much larger systems. For example, when A = 1000 the relaxation times estimated by 
Henon's value and the value 7 = 0.02 differ by about 25%. Even the largest of our models 
are, therefore, a somewhat uncertain guide to the absolute time-dependence of relaxation 
processes in larger iV-body systems. 

The post-collapse evolution is, to a first approximation, self-similar, as one can see 
from the fact that the curves in Fig.l are very nearly parallel, in the sense that all can 
be brought nearly into coincidence by a suitable vertical displacement. This is evidently 
not true in the early post-collapse phase for the outermost radii (mass fractions of 75% or 
more). These aspects of the evolution are illustrated qualitatively in Fig. 2, which shows 
the surface density profile scaled (radially) to the half-mass radius. In post-collapse the 
profiles could be brought into near-coincidence by a vertical shift. Nevertheless even for 
the latest profiles shown in Fig. 2 the outermost radius at which the density is plotted 
(which corresponds to a mass fraction of 90%) moves outwards slowly as time advances, 
indicating a modest but sustained departure from homologous evolution. 

Though there are departures from self-similarity, the half-mass radius still varies in 
approximate accord with naive theoretical expectations, i.e. oc t 2 / 3 . In fact there are 
small departures from this formula. To clarify this point we carried out a similar exercise 
to that described in Paper II, i.e. we fitted the following functions to the A-body data: 

N = a(t-to)- v , 

r h ^b(t-toY 2+ ^ 3 , (1) 
(Nrl/m) 1 ' 2 = c(t-t Q ), 

where A is the number of bound stars, is the half-mass radius, t is time, m is the mean 
individual mass of stars inside r^, and a, 6, c, v and to are fitting parameters. The values 
obtained are given in Table 3. 

According to standard relaxation theory (Spitzer 1987) the coefficient c should be 
proportional to the Coulomb Logarithm, usually taken as ln(7iV) for some value of the 
constant 7. But as can be seen from results presented in Table 3 it is a good approximation 
to assume that c is constant, i.e. independent of A, and this is inconsistent with any 
plausible value of 7. A similar conclusion follows from consideration of the values of b. We 
have no definitive explanation for these results, though it should not be concluded that 
the evolution of systems of stars with unequal masses behaves in a manner inconsistent 
with standard relaxation theory. At some level of detail the mass distribution evolves 
differently in systems with different A, and also our results have been obtained by fitting 
to the evolution over an A-dependent number of relaxation times; for example, the factor 
by which expanded was very different in the different models. These A-dependent effects 
may have conspired to mask the expected A-dependence of the Coulomb logarithm. 

It is to be noted that Fig.l illustrates projected Lagrangian radii, i.e. containing a 
given fraction of the projected mass. The corresponding picture for the conventional (three- 
dimensional) Lagrangian radii is qualitatively similar, at least in respect of the aspects to 
which attention has been drawn above. The three-dimensional half-mass radius is larger in 
the post-collapse phase by about 30%. In the same way the space density profiles exhibit 
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very similar features to the surface density profiles shown in Fig. 2. Also the results for 
N = 1000 and 250 do not differ qualitatively from those shown here for N = 500. The 
main quantitative difference, which is the time scale of the evolution, has already been 
described above. 

3.2 Mass segregation 

At first sight it is surprising that the structure of the cluster in post-collapse evolution 
appears to be so nearly self-similar, and that the naive theoretical expectation, eq.(l), is 
so well satisfied, because there is no reason for supposing that mass segregation should 
have a negligible effect after core collapse. Therefore, even though it remains true that 
the time scale for post collapse expansion should be governed by a time of relaxation, it is 
not as clear as it is in the case of equal masses that this time scale may be characterised 
by the half-mass relaxation time: as conventionally defined (Spitzer 1987) this depends 
on the mean mass as well as the total mass and the half-mass radius, and one might 
have expected the somewhat different time scale for mass segregation to complicate the 
evolution of the half-mass radius. In the present section we shall discuss mass segregation, 
and we shall find an empirical explanation for the fact that the evolution of the half-mass 
radius is so simple. What is still lacking, however, is a theoretical understanding of why 
mass segregation in post-collapse evolution takes the simple form observed. 

The basic result is illustrated in Fig. 3. Extremely rapid segregation of masses takes 
place throughout core collapse. The remarkable conclusion to be drawn from this figure, 
however, is that the profile of mean masses evolves remarkably little throughout the re- 
mainder of the calculations, which amounts altogether to about 20 collapse times. It is 
presumably this fact, i.e. the fact that the distribution of stellar masses is nearly static 
(in Lagrangian coordinates), which helps to explain why mass segregation does not, as one 
might expect, complicate the post-collapse evolution of the cluster. Note, however, that 
this merely places the problem one step back: we have no explanation for the fact that 
the distribution of masses at each Lagrangian radius shows such little evolution after the 
initial phase of mass segregation. 

While Fig. 3 illustrates the projected distribution of the mean mass, the spatial distri- 
bution shows very similar features. There is, however, a slight but noticeable drop in the 
mean mass inside the 20% radius in the early post-collapse regime and for all Lagrangian 
radii during the advanced post-collapse phase. This drop is more pronounced than the 
corresponding drop in the projected data (Fig. 3), but still no greater than about 0.1 dex. 
It is not clear that this is correctly attributed to mass segregation, since stars of high mass 
are removed from the core and the whole system by binary formation, evolution and escape 
(§6). 

Since the distribution of masses is almost constant in the post-collapse regime, except 
for the overall expansion of the radii and the corresponding decrease in the density, we 
illustrate the density distribution of the different masses, binned as stated in §2, at one 
point in the post-collapse phase (Fig. 4). Initially, of course, all density profiles would 
be identical, except for an overall vertical displacement from one species to another. As 
can be expected at this point in the post-collapse evolution, however, the effects of mass 
segregation are clearly visible. The most massive stars dominate the central parts of the 
systems and the less massive stars dominate the halo. As a rule, the larger the stellar 
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mass the greater the concentration to the centre, with one exception: stars from the most 
massive bin are not, as one might have expected, confined to the central and middle parts of 
the system; they are present in the outer halo as well, at distances even greater than those 
of the least massive stars. (More precisely, their 90% Lagrangian radius exceeds that of 
the least massive stars.) This behaviour is probably connected with the ejection of massive 
stars and binaries from the core to the outer parts of the system during interactions between 
binaries and field stars. This is evidence that binary activity can disturb the otherwise 
clear picture of mass segregation in these systems. 

4 EVOLUTION OF THE VELOCITY DISTRIBUTION 

4.1 The overall velocity dispersion 

Many aspects of the overall evolution of the velocity distribution are reflected in Fig. 5. 
After the relatively brief period of core collapse the decline in the mass-weighted root 
mean square projected speed is nearly homologous. The profiles of the projected velocity 
dispersion behave in a manner very similar to that of the density profiles (cf. Fig. 2): if the 
radial coordinate is scaled by the projected half-mass radius the profiles at different times 
differ from each other only by a vertical shift in a logarithmic plot, to good approximation. 
This is also true of the three-dimensional velocity dispersion profiles, though they are very 
slightly steeper: between the centre and the 90% Lagrangian radius the difference in rms 
speed is about 0.36 dex for the three-dimensional profile, and about 0.33 for the projected 
profile. 

4.2 Equipartition 

The most interesting kinematic distinction between the models discussed here and 
those analyzed in Papers I and II concerns the question of equipartition, or more generally 
how the velocity dispersion differs from one population to another. Just as we have at- 
tempted to summarise the distribution of masses by the mean mass profile (Fig. 3) we may 
begin by attempting to define a single measure of the extent to which different masses are 
in equipartition, as follows. 

For a given Lagrangian shell let Tj~ be the mean kinetic energy of the stars in that 
shell in mass class k, with 1 < k < 10 (§2), and T the mean kinetic energy of all stars in 
the shell. If there happen to be no stars in a given mass class in that shell then we define 
7^ = 0. Next we define the 'equipartition parameter' as 

B = (2) 

where iV> is the number of mass classes for which Tk > 0. The obvious feature of B 
is that in the case of equipartition (of kinetic energies) we should have = T for all 
k, and so B = 1. Likewise, in the case of equipartition of velocities we would have 
B = ^mfe/(10(m)), where nik is the mean mass of stars in the kth. mass class. Notice 
that X] m fc/(10( m )) 7^ 1, since the average on the left hand side is not weighted by the 
number of stars in each mass class. For our initial mass function, for example, we would 
have B = 3.8 (i.e. logS = 0.58) in the case of equipartition of velocities. 
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Clearly one may compute B either for three-dimensional Lagrangian shells or for 
projected annuli, in which case could be defined as the mean kinetic energy of motion 
in the line of sight. This has been adopted for Fig. 6, which shows profiles at various times. 
The curves are noisier than those previously shown, because the small number of stars of 
high mass have relatively high weight in the numerator. Nevertheless one can easily observe 
an extremely rapid evolution (in the first 5 time units for N = 500) from the initial flat 
profile to one which remains almost constant throughout the remainder of core collapse 
and in the post-collapse phase provided (as has been done here) the radial coordinate is 
scaled by the projected half-mass radius. In the core the value is close to that expected for 
equipartition of energy, whereas in the outermost parts of the system the departure from 
energy equipartition, to the extent that this is measured by £?, is even greater than it was 
initially. 

More detail is shown in Fig. 7, which gives the profiles of the mean square projected 
velocity as a function of projected radius, for all mass groups. (We give this for only 
one time, since Fig. 6 suggests that there is little evolution after an initial short period 
of adjustment.) Unfortunately, because of the small number of stars in each mass bin 
(particularly the massive bins) the data presented in Fig. 7 are very noisy. However, the 
main features of this picture are clearly visible. In the core, dominated by massive stars, 
the mean square projected velocities follow closely the rule of energy equipartition. In the 
outer halo, dominated by less massive stars, the mean square projected velocities behave 
differently: they are nearly the same for most of the mass bins. This reflects the initial 
conditions, for which all mass bins were in velocity equipartition. In the outer halo the 
relaxation time is very long (compared to that in the core, or even at the half-mass radius), 
and so the evolution is very slow and the initial conditions are well preserved. Therefore 
the equipartition parameter B increases towards the outer parts of the system (see Fig. 6). 

For the most massive bin the velocity is nearly constant throughout the system. The 
explanation for this may be connected again with binary activity, which removes stars and 
binaries from the core and ejects them into the halo with velocities high compared with 
most other stars in the halo. Because massive stars are preferentially removed from the 
core the behaviour of the most massive bin is markedly different from that of the others. 

4.3 Anisotropy 

In discussing the anisotropy we shall concentrate on projected values (Fig. 8). Thus 
each annulus is bounded by circles containing fixed fractions of the total mass (relative to 
the density centre of each model), and within each annulus we denote by (v%) and (v^), 
respectively, the mean square radial and transverse components of the projected velocity. 
Then the anisotropy, defined by (v^)/(v^), is plotted at a radius given by the mean of the 
radii of the bounding circles scaled by the projected half-mass radius. 

It can be seen from Fig. 8 that the evolution of the anisotropy profile occurs on a 
longer timescale than that of mass segregation (Fig. 3) and the tendency to equipartition 
(Fig. 6). (This difference of timescales is not really surprising, since the main evolution of 
the anisotropy occurs at large radii, unlike the other two phenomena.) When N = 500 
the profile continues to evolve in the outer parts until about t ~ 100, and only thereafter 
remains relatively constant (when the radius is scaled by the expanding projected half-mass 
radius). 



8 



As already stated the values plotted in Fig. 8 are projected values. In the quasi-steady 
phase after t ~ 100 the three-dimensional values of log((vf)/(vf)) are about —0.6 at the 
outermost radii shown on Fig. 8 (i.e. at a similar value of the three-dimensional radius 
scaled by the three-dimensional half-mass radius). Though this might be taken to imply 
that the three-dimensional anisotropy is smaller, it must be recalled that the value of the 
logarithm in an isotropic distribution is log 2 in three dimensions, and for projected data. 

Another point to be made about these data concerns the iV-dependence of the aniso- 
tropy. As in the case of equal masses (Giersz & Spurzem 1994) the anisotropy in the outer 
parts saturates at a value which is greater (i.e. more radially anisotropic) for smaller N. 
For N = 250, 500 and 1000 the values of the logarithm (three-dimensional anisotropy) are 
about —0.85, —0.63 and —0.48, respectively. This may have a similar explanation to that 
of the faster expansion of the outer Lagrangian radii for smaller N in systems with equal 
masses (cf. Paper II), i.e. the fact that the relatively shallower potential well of systems 
with small N makes it easier for stars which are ejected from the core to populate the far 
halo. 

The way in which the anisotropy is distributed among the different mass classes is 
illustrated in Fig. 9, which exhibits the projected anisotropy profile for N = 500 at t = 300, 
i.e. well into the phase in which the overall profile appears to have settled down into 
equilibrium (as a function of radius scaled by the half- mass radius). Generally, the core 
is isotropic for all mass bins (with very substantial fluctuations for massive bins). In the 
outer halo the anisotropy is bigger for the bins corresponding to the most massive stars; as 
with equipartition and the velocity dispersion, the most massive bin is exceptional. This is 
again consistent with the suggestion that the massive stars are ejected from the core into 
the halo by binary activity. 

5 EVOLUTION OF THE CORE 

In systems of stars of equal mass there are already two commonly used measures of 
the core radius. One, stemming from the work of King (1966), defines the core radius r c 
in terms of the central one-dimensional velocity dispersion a c and the central mass-density 
p c by the formula 

< 3 » 

The other, derived from the work of Casertano & Hut (1985), defines it as 

rl = ^M, (4) 



where pi is a measure of the density at the location of the ith star, and is its distance 
from the density centre (§2). (The formula proposed by Casertano & Hut actually differed 
from eq.(4) in that all the terms were unsquared.) In our work we have estimated a c 
and p c from the innermost 1% by mass (cf. below), while pi was estimated as in §2. It 
is known (cf. Aarseth & Heggie 1992) that the definition in eq.(3) gives a larger value 
than eq.(4) when applied to a Plummer model, which is the starting configuration of our 
calculations. But in our unequal-mass models it quickly becomes the smaller, exhibiting a 
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much deeper collapse: for N = 500, for instance, the core radius defined by eq.(3) decreases 
by a factor of 3 approximately, whereas the decrease for r c defined by eq.(4) is only by 
a factor of about 1.4. Throughout the post-collapse evolution eq.(4) gives a core radius 
which is larger than that defined by eq.(3) by a fairly consistent factor which is nearly 2. 
We prefer eq.(3), entirely on the grounds that it has a dynamical significance, which has 
not been demonstrated for eq.(4). 

As for the value of r c itself, the simplest way of expressing the result for post-collapse is 
to give the ratio to the half-mass radius (Table 2), which is almost constant in this phase 
of the evolution. These results are quite interesting, as they contradict the theoretical 
expectation for equal masses, which implies that, for large N, the ratio should vary as 
iV -2 / 3 (Goodman 1987). Whether this implies that the models have not yet reached the 
asymptotic behaviour of systems of sufficiently large N, that systems with stars of unequal 
mass behave in a qualitatively different way, or that our adopted definition of core radius 
is faulty in some way for systems with unequal masses, cannot be decided on the basis of 
our data. Certainly, our conclusion for systems of equal mass (Paper II) was that the size 
of the core was roughly in line with theoretical expectations, even for the small values of 
N that were studied. 

The definition of core radius having been decided on, we can readily determine the 
number of stars within this radius (Table 2), from which it appears that N c ~ N/60 
throughout almost all of the post-collapse evolution. The mass of stars within r c is ap- 
proximately 0.05, independently of N, and so the mean stellar mass in the core is 3/iV 
approximately, i.e. about three times the initial global mean mass (cf. also Fig. 3, which 
gives a similar value for the projected central mean mass). 

Our estimate for the central velocity dispersion, which we now discuss, is obtained 
from the mass- weighted mean square speed of stars within the innermost 1% by mass. Even 
though we include all stars from all models (at given N), because of mass segregation the 
stars in this zone are relatively massive and few in number. Therefore the estimate of 
the mean square speed is subject to quite large fluctuations, especially in relation to the 
rather small changes in this quantity during the evolution. Fig. 10 shows the results for 
all values of N that we have studied. One reason why these results are interesting is the 
initial decline in the early phase of core collapse, though because of the fluctuations the 
result is convincing really only for N = 1000. A similar drop was noticed in the equal-mass 
calculations discussed in Paper I, but it is more natural here to interpret it as a consequence 
of mass segregation and the tendency towards equipartition. A similar drop can be seen 
as well in data presented by Chernoff & Weinberg (1990) for Fokker-Planck calculations 
for King models with unequal masses. As already stated, and discussed further below, the 
core quickly becomes dominated by massive stars, whose rms speed drops below that of 
the stars of lower mass. Thereafter the core collapse, which involves mainly the dominant 
massive stars, leads to the expected increase in the rms speed. 

The fact that the central density rapidly becomes dominated by the most massive 
stars is responsible for the evolution of the mean mass within the innermost Lagrangian 
radius, as shown in Fig. 3. Fig. 11 gives more detail, showing the contribution to the central 
projected density from the several mass bins. (Note, however, that the density is influenced 
also by the sizes of the mass groups, which are uniform in logra.) 
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Another core parameter of importance (e.g. for escape of binaries, and the efficiency of 
heating by three-body binaries) is the depth of the potential well at the centre of the system 
(Fig. 12). The minimum is less sharp than in the case of equal masses (Paper II), perhaps 
because of an even greater relative spread in the times of core collapse in the models with 
unequal masses, and perhaps because the presence of a few massive stars in the core causes 
larger fluctuations than in the case of equal masses. Note that the minimum depth of the 
potential is increasing with N. This result is quite interesting as it contradicts theoretical 
expectation (e.g. Goodman 1987) and numerical results for equal masses (Paper II), which 
implies that the minimum of the potential decreases with increasing N. We interpret this 
as a richness effect: the maximum individual stellar mass which is present in a given small 
iV-body system may be much less than the maximum allowed mass, if the mass spectrum 
is steep. (There is evidence for this in Table 1, which shows that the fraction of mass in 
the heaviest mass group is an increasing function of N.) One consequence of this is that 
a small iV-body system is unlikely to be able to form a binary from stars close to the 
maximum allowed mass, and therefore binary heating is likely to be less efficient than in 
a system with sufficiently large N. 

All the results on core evolution displayed so far concern the time dependence, but 
it is also instructive to display the evolution of two fundamental core parameters against 
each other. For Fig. 13 we have chosen the central mass-density and the central mass- 
weighted root mean square speed. This bears some resemblance to the corresponding 
result for equal masses (Paper II). The dip in the central velocity dispersion during core 
collapse has already been discussed, but the post-collapse behaviour is another example 
of the surprisingly simple behaviour of these unequal-mass models in this phase of the 
evolution: if it is assumed that p c and v\ are proportional to values at the half-mass 
radius, then it follows that p c i>~ 6 oc M~ 2 . Since the total mass M varies slowly (cf. §7), 
this predicts an approximate power-law relation which is quite closely verified in Fig. 13. 
It is quite surprising that the presence of a spectrum of masses does not invalidate such 
simple arguments. 

6 BINARIES 
6.1 Definitions 

For systems of stars with equal masses the notions of "hard" and "soft" binaries 
are relatively uncomplicated. When there is a spectrum of masses, however, it becomes 
less clear what definition to adopt. The simplest approach would be to classify binaries 
according to the ratio ej (kT), where e is the (internal) binding of the binary, and 1.5/cT is 
the mean kinetic energy of all stars in the neighbourhood of the binary. On the other hand 
it is not clear that this is the most useful definition when the masses may be very unequal. 
In the case of equal masses the notions of "soft" and "hard" are important usually because 
they correspond, roughly, to the distinction between binaries which tend to be destroyed 
by encounters and those which become harder (Gurevich & Levin 1950). While work by 
Hills (1990) and by Sigurdsson & Phinney (1993) gives much useful information on the 
effect of encounters between a binary and a single star when the masses are different, it is 
still difficult to decide whether a binary, with a given energy and component masses, will 
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tend to harden or disrupt under the action of encounters with a field of single stars, when 
the spectrum of stellar masses may be very time- and position-dependent. 

In view of this difficulty we have adopted a definition which is much simpler to imple- 
ment, but is more difficult to interpret physically. In what follows we define a "hard" binary 
to be a regularised binary, in the sense of the code NB0DY5 (Aarseth 1985). We shall 
also refer more briefly to "soft" binaries, which are defined to be unregularised binaries 
with binding energies greater than 0.1/cT (cf. §2). 

Clearly, many binaries fall into neither of the classes which we term hard and soft, 
because there will in general be a number of binaries whose internal binding energies are 
too low to satisfy the criteria for regularisation. There is no evidence, however, for the 
presence of long-lived binaries which do not fall into the class of regularised pairs. 

6.2 Statistics of Binaries 

As with equal masses (Paper II) the onset of the formation of hard binaries is quite 
abrupt (Fig. 14). Though this is to be expected, being associated with the end of core 
collapse, what is perhaps more surprising is the result in the post-collapse phase, where the 
number of hard binaries is nearly constant in time, and independent of A. In this phase the 
mean value is about 1.3. For systems of equal masses (Paper II) the corresponding result 
is about 2.4 (but there are some indications for A-dependence of the form: Nb oc A 1 / 5 6 ). 

The spatial distribution of these binaries is heavily centrally concentrated: in the early 
post-collapse phase about 80% of the hard binaries lie within the radius containing 20% 
of the mass of all stars. The central concentration of hard binaries decreases somewhat in 
later phases of post-collapse expansion, reaching a nearly constant value (as measured by 
the number of binaries within the 20% Lagrangian radius) for times greater than about 
100 A-body time units (which apparently coincides roughly with the time at which the 
anisotropy levels off). This drop is accompanied by an increase in the number of hard 
binaries in the outer parts of the cluster. This is particularly clearly visible in data for 
A = 250, which extends to a larger number of collapse times than our data for the larger 
systems. 

The total internal energy of all bound hard pairs also shows the expected sharp change 
towards the close of core collapse (Fig. 15), and in the early post-collapse phase it follows 
closely a power-law dependence, which can be interpreted very roughly as follows. It is 
known from the behaviour of the half-mass radius (§3.1 above) that the "external" energy 
of the bound stars (i.e. excluding the internal energy of bound binaries) varies roughly as 
E b ext ~ —C(t — to) _2 ^ 3 , where C is a positive constant, though this form of dependence 
would be modified slightly if account were taken of mass-loss. If no binaries escape and 
we neglect the flux of energy associated with escapers, it follows that the internal energy 
of bound binaries varies roughly as 

El t ^E + C(t-t )- 2 /\ (5) 

where Eq is the initial energy of the system, i.e. —0.25 in our units. 

In fact eq.(5) provides a reasonable fit to the data early in the post-collapse phase, 
until one of the above assumptions begins to break down seriously. Indeed the rise in E h int 
after this early phase indicates the escape of hard binaries (cf.§7), and then the derivation 
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of eq.(5) breaks down. Before this happens, however, we see from Fig. 15 that E b int falls 
below the minimum value predicted by eq.(5), and the reason for this is that escape by 
two-body encounters (§7) tends to decrease the external energy of the system, i.e. to make 
it more bound, thus decreasing the effective value of Eq. 

The mean mass of the hard binaries is also instructive (Fig. 16). In the innermost parts 
of the cluster they tend to consist of the most massive stars, whose mean mass declines 
slowly with time as the most energetic (and most massive) binaries escape or are ejected to 
the outer parts of the cluster (§7). Indeed this gives rise to the perhaps rather surprising 
behaviour of the hard binaries in the outer parts, whose mean mass tends, if anything, to 
increase. It is in the intermediate radii that the mean masses of binaries are smaller, and 
indeed decline relatively rapidly as the evolution proceeds. 

The assertion that the massive binaries in the outer parts are ejected from the core 
gains further support from an examination of the eccentricities of the hard binaries. Out- 
side the half-mass radius the mean eccentricity of hard binaries, which is about 0.8, sig- 
nificantly exceeds that of the hard binaries inside (about 0.7, i.e. close to the thermal 
average value of 2/3 (Jeans 1929)). 

7 GLOBAL EVOLUTION 
7.1 Escape 

In the previous section attention was already drawn to the behaviour of massive hard 
binaries, which tend to be ejected into the outer parts of the cluster if not out of the cluster 
altogether. Surprisingly, perhaps, this escape mechanism is sufficiently effective that the 
mean mass of the remaining bound members of the cluster tends to decrease, whereas one 
might have thought that the preferential escape of low-mass stars in two-body encounters 
(Henon 1969) would lead to an increase in the mean mass. This is illustrated in Fig. 17, 
which shows that the mass of the cluster decreases by 15% at a time when the number of 
stars has decreased by less than 11%. 

These mean data actually mask a more complicated picture, which is given in Fig. 18. 
At early times, i.e. during core collapse (Fig. 18a) there is indeed preferential escape of 
the stars of lowest mass, as might be expected from Henon's theory and as a byproduct of 
mass segregation. The actual fractional escape is less than 1% however. After core collapse 
(Fig. 18b) the mass-dependence is reversed, as the most massive stars preferentially escape. 
As we have seen (Fig. 3) there is little further mass segregation after core collapse, and 
therefore part of the mechanism which led to the preferential escape of the lightest stars is 
suppressed. Now, however, interactions in the core, which consists predominantly of stars 
of high mass, some in binaries, lead to escape of stars directly from the core, and this 
depletes the population of stars of highest mass. 

The number of binary escapers is very small. For example when N = 500 the average 
number is about 0.5 at t ~ 400, ~ 20t co u (cf. Table 2); this number is about 1% of 
all escapers (cf. Fig. 17). Simple theory implies that the possibility of escape of binaries 
depends, among other things, on the mass of stars which form the binary. The fraction 
of the kinetic energy which the centre of mass of the binary (of total mass m^) can gain 
during an interaction with a field star of mass ms is approximately equal to m^/ {mb + m^). 
Therefore when m& > 2m^ this fraction is smaller than in the case of equal masses. It 



13 



follows that a massive binary has to experience more (or more energetic) interactions with 
field stars before its removal from the system. In a sense, therefore, it remains in the core 
longer than a binary which consists of stars with masses equal to that of the single stars. 

7.2 The energy budget 

Initially the energy of the iV-body systems consists entirely of the "external" energy 
of the bound members, which we denote by E ext . When binaries form, E ext increases 
(i.e. the cluster becomes less bound), and the "internal" energy of the bound binaries 
E h int decreases in compensation (Fig. 15). At the same time the escape of single stars by 
two-body processes diverts energy into the external energy of the escapers E^ xt , and this 
tends to decrease the energy E ext of the remaining bound stars. Finally, as interactions 
with binaries become energetic enough to lead to their escape (the energy associated with 
the motion of their barycentres being denoted by E™. t ), the internal energy of escaping 
binaries (Ef% t ) decreases (i.e. increases in magnitude), and the internal energy of bound 
binaries E h int increases (decreases in magnitude) in compensation. The balance between 
these various processes is shown in Fig. 19. Note that, although the number of binary 
escapers is small (§7.1) they have a very significant effect on the energy budget. It is 
worth comparing Fig. 19 with the corresponding diagrams published by Heggie & Aarseth 
(1992) for larger systems (N ~ 2500) containing primordial binaries, though in Fig. 19 we 
have split E^ xt into two components, corresponding, respectively, to the kinetic energy 
of single stars and binaries. Though the initial value of E int is non-zero in that case, 
the trends in the data are remarkably similar, even though the calculations of Heggie & 
Aarseth were larger and used equal masses. 

8 DISCUSSION AND CONCLUSIONS 

We have analysed the evolution of isolated iV-body systems in which the spectrum 
of stellar masses is a power law. The systems contained between 250 and 1000 stars, and 
for each iV we have averaged results from a large number of runs, in order to improve the 
statistical quality of the data. 

Our results confirm the extreme rapidity of core collapse in the presence of a mass 
spectrum (cf. Inagaki & Wiyanto 1984). Theoretically, the iV-dependence of the rate of 
core collapse involves the Coulomb logarithm In "fN, but our results do not constrain the 
value of 7 as strictly as in the case of equal masses (Paper I). Our results are consistent with 
a value around 7 ~ 0.02, which is much smaller than the corresponding result for equal 
masses, in qualitative agreement with the theory of Henon (1975). For such relatively small 
values of N as in our simulations the interpretation of the time scales is quite sensitive to 
the value of 7. For example the value of the initial half- mass relaxation time t r h(0) for 
N = 500 varies from about 20 for 7 = 0.02 to 11 for the standard value 7 = 0.4. If the 
former value is adopted we find that core collapse is complete (Table 2) in about 0.9^(0). 
Lest this seem much too short, it may be noted that Chernoff & Weinberg (1990), using 
a Fokker-Planck code, found that the time to core collapse was also about 0.9t r / l (0) for 
the same mass spectrum, though the initial model they adopted was a King model with 
dimensionless central potential Wo = 3, and a tidal cutoff was included. 

One of our most curious findings (§3.2 and Fig. 3) is that the spatial distribution of the 
mean stellar mass is nearly constant (i.e. time- independent) throughout the post-collapse 



14 



evolution, when the overall expansion of the system is scaled out. While this helps to 
explain why the time-dependence of the post-collapse expansion takes the simple form of 
eq.(l), it is only an empirical explanation. It remains to be understood on theoretical 
grounds why the extremely rapid segregation of masses during core collapse comes to such 
an abrupt halt. One way in which this question might be approached would be to try to 
find a self-similar, multi-mass, post-collapse solution of the Fokker-Planck or gas equations, 
in analogy with the equal- mass models of Goodman (1984, 1987) or Inagaki & Lynden- 
Bell (1983). It was suggested long ago by Lynden-Bell (pers. comm.) that one might 
search for such a solution for core collapse, in which the distribution of stellar masses was 
a power law. Be that as it may, there is now empirical evidence that a search for such a 
post-collapse solution might be fruitful. 

It is equally surprising, perhaps, that the departures from equipartition (Fig. 6) also 
follow a very simple behaviour in the post-collapse regime. Certainly it may be expected 
that departures from energy equipartition will be small in the core, where the relaxation 
time is short, and this is what is found. But it might have been expected that, at radii 
where the departure from equipartition becomes significant, one should have observed a 
progressive tendency towards equipartition. However, we have found that there is essen- 
tially no change in the profile of our equipartition parameter, except for that caused by 
the nearly homologous expansion of the entire system. 

Another surprise, except in hindsight, is the evolution of the heaviest stars. We have 
seen (§7.1) that they escape at a relatively high rate, which would not necessarily be 
expected because it might be thought that they would tend to eject stars of low mass 
while settling, through mass segregation, into the core. Perhaps this does happen, but in 
the core the massive stars tend to form energetic binaries, whose further interactions tend 
to eject these binaries and the moderately massive stars in the core, thus enhancing their 
rate of escape. 

Some hard binaries are ejected into the outer parts of the cluster without completely 
escaping, and we noticed (§6.2) that these tend to have higher eccentricities than average, 
presumably as a result of the ejecting encounter. Thus hard binaries beyond the half-mass 
radius tend to have large masses and large eccentricities. Though our models obviously 
differ in crucial respects from the galactic globular clusters, it is possible to see here 
an example of the dynamical behaviour which gives rise to such objects as the binary 
millisecond pulsar in M15 (Phinney 1993). 
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Table 1 

MASS GROUPS 



range of 


masses 


mean mass 






(unit 


s of TTlrain) 












N= 


250 


1.000 


1.439 


1.186 


0.1895 


101.1 


1.439 


2.067 


1.706 


0.1722 


63.63 


2.067 


2.970 


2.454 


0.1428 


36.71 


2.970 


4.267 


3.527 


0.1161 


20.84 


4.267 


6.132 


5.054 


0.09463 


11.85 


6.132 


8.810 


7.213 


0.08564 


7.446 


8.810 


12.66 


10.34 


0.06648 


4.054 


12.66 


18.19 


15.01 


0.06046 


2.554 


18.19 


26.13 


21.36 


0.04526 


1.321 


26.13 


37.50 


31.18 


0.02588 


0.5357 



initial data 
mass fraction and number 



IN = 


:O00 


1\T 

IN = 


1000 


0.1902 


202.4 


0.1987 


422.9 


0.1705 


126.1 


0.1721 


254.3 


0.1421 


73.07 


0.1440 


148.5 


0.1198 


42.88 


0.1210 


86.55 


0.09874 


24.66 


0.09822 


49.03 


0.07981 


13.96 


0.08138 


28.35 


0.06978 


8.518 


0.07357 


17.85 


0.05607 


4.714 


0.05732 


9.675 


0.04412 


2.607 


0.04220 


4.975 


0.02777 


1.125 


0.03547 


2.900 



Note: the mass fractions and numbers were determined from the sets of initial conditions. 
For an analytical mass function, the mass fraction would be independent of N, and the 
number would be proportional to N. 
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Table 2 



SUMMARY OF RESULTS 

N t coU r c /r h N c 

250 12 (95) 0.11 4 

500 18 (164) 0.14 8 

1000 21 (332) 0.16 17 



Note: tcoii is the approximate time of the end of core collapse (determined by the minimum 
of the central potential, cf. fig. 12). Values for equal masses (from the data discussed in 
Papers I and II) are shown in brackets. The values of r c /rh and N c are representative of 
almost all of the post-collapse evolution. 

Table 3 

FITTED PARAMETERS 



N t a 

250 -5.688 316.94 

500 -8.256 622.95 

1000 +4.661 1236.6 



v b c 

-0.067 0.159 16.95 

-0.063 0.095 17.24 

-0.061 0.064 17.25 



Note: the symbols are explained in eqs.(l) and the accompanying discussion. 



18 



FIGURE CAPTIONS 



Fig.l Evolution of projected Lagrangian radii for iV = 500. Also shown is the (three- 
dimensional) core radius, defined in §5. 

Fig. 2 Evolution of the projected (surface) mass density profile for A = 500. The radial 
coordinate is scaled by the instantaneous projected half-mass radius. The densities plotted 
are average values between each Lagrangian radius, and have been plotted at the mean 
radius in each Lagrangian shell. 

Fig. 3 Evolution of the mean individual stellar mass within each projected Lagrangian shell 
for A = 500. Thus the curve labeled "2%" gives the mean masses of stars whose projected 
distance from the density centre lies between projected radii which enclose 1% and 2% of 
the projected mass, respectively. A binary was treated as a single star with a mass equal 
to the total mass of the components. 

Fig.4 Profile of projected mass density of each mass class at time t = 300 for A = 500. 
The Lagrangian radii have been determined separately for each mass class, and since the 
plotted values of £ represent average estimates for annuli separated by these projected 
Lagrangian radii, the innermost and outermost plotted points vary from one species to 
another. 

Fig. 5 Root mean square projected velocity in each projected Lagrangian annulus for A = 
1000. The data have been cleaned by removal of energetic escaping particles, which would 
otherwise cause a "spike" , especially in the outer shells, where they spend the longest time. 
Also plotted (v c ) is the rms projected velocity within a sphere whose radius equals the core 
radius. (Note that this quantity is not a projected velocity dispersion in the usual sense.) 
Fig. 6 Projected equipartition parameter B (eq.(2)) at various times in a set of models 
with A = 500. The radial coordinate has been scaled by the projected half-mass radius. 
Fig. 7 Projected velocity dispersion profiles for each mass group for A = 500 at time 
t = 300. 

Fig. 8 Projected anisotropy profiles for A = 500. 

Fig.9 Projected anisotropy profiles at t = 300 for each mass class in models with A = 500. 
Fig. 10 Mass- weighted root mean square central speed as a function of time. The result 
for A = 250 is plotted to scale, while those for A = 500 and 1000 have been displaced 
vertically by 0.3 and 0.6 units, respectively. The data have been smoothed over two time 
units. 

Fig. 11 The projected mass density in the eight mass groups within the innermost La- 
grangian shell (i.e. the 1% Lagrangian radius), for A = 500. 

Fig. 12 Evolution of the central potential (i.e. at the location of the density centre) for 
A = 250, 500 and 1000. 

Fig. 13 Evolution of core parameters p c and v c for A = 500. The straight line corresponds 
to p c oc v®- The A-body data were smoothed using the standard routine SMOOFT from 
Press et al (1986). 

Fig. 14 Number of hard binaries as a function of time, for A = 250, 500 and 1000. Only 
binaries which are not regarded as escapers are included. 

Fig. 15 Total internal energy of all bound hard binaries for A = 250, 500 and 1000. 
Fig. 16 Mean mass of hard binaries in 4 spatial zones for A = 500. M to t is the total bound 
mass of the system. 
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Fig. 17 Evolution of total number and mass of bound stars in systems with N = 500 stars 
initially, expressed fraction of the initial value. 

Fig. 18 Fractional escape of stars in different mass classes, for N = 1000. Fig. 18a shows 
details for the collapse phase and the start of the post-collapse evolution, while Fig. 18b 
shows the long-term trend. 

Fig. 19 Energy balance for N = 1000. The meaning of the symbols is given in the text. 
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